home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / demos / npend / maple / euler.map < prev    next >
Text File  |  1999-09-16  |  4KB  |  165 lines
  1. #----------------------------------------------------------------------------
  2. #     GENERAL PART
  3. #----------------------------------------------------------------------------
  4. #----------------------------------------------------------------------------
  5. # Notation for Tex output 
  6. # and TeX output routines
  7. # this routines generates proper derivative notations up to second derivative
  8. # for a set of variable. The TeX name for the variables is also given
  9. # Warning this function will detroy previous values of texsub 
  10. # Ex : addnotations( { x = `x`  , th = `\\theta`, y =`y` } );
  11. #-----------------------------------------------------------------------------
  12.  
  13. addnotations:=proc(varlist)
  14.     local x,nots:
  15.     nots:={};
  16.     for x in varlist do
  17.         nots:={op(nots),x, 
  18.             cat(lhs(x),`d`)=cat(`\\dot{`,rhs(x),`}`),
  19.             cat(lhs(x),`dd`)=cat(`\\ddot{`,rhs(x),`}`)}
  20.     od;
  21.     texsub:=nots
  22. end:
  23.  
  24. sortieinit:=proc(filename)
  25.     writeto(filename);
  26.     writeto(terminal);
  27. end:
  28.  
  29. sorties:=proc(filename,comment,exp)
  30.     appendto(filename);
  31.     lprint(comment);
  32.     lprint(` `);
  33.     mtex(exp);
  34.     writeto(terminal);
  35. end:
  36.  
  37. sortiesI:=proc(filename,comment)
  38.     appendto(filename);
  39.     lprint(comment);
  40.     writeto(terminal);
  41. end:
  42.  
  43. sortiesM:=proc(filename,comment,expr)
  44.     appendto(filename);
  45.     lprint(comment);
  46.     lprint(` `);
  47.     map(x -> mtex(x),expr);
  48.     writeto(terminal);
  49. end:
  50.  
  51. readlib(tex):
  52. texwid:=120;
  53. sortieinit(`systeme.tex`);
  54.  
  55. #---------------------------------------------------------------
  56. # functions for computing euler equations 
  57. # L : the Lagrangian,
  58. # q,qd,qdd are the state vector and its derivatives
  59. #---------------------------------------------------------------
  60. with(`linalg`):
  61.  
  62. euler_equations:=proc(L,q,qd,qdd)
  63.      local k,m:
  64.     m:=nops(q);
  65.     v:=matrix(m,1,0);
  66.     for i to m  do  
  67.            v[i,1]:=LL(q[i])=simplify(time_diff(diff(L,qd[i]),q,qd,qdd)
  68.                 -diff(L,q[i]));
  69.     od;
  70.     eval(v):
  71.     end:
  72.  
  73. #---------------------------------------------------------------
  74. # Time derivative computation of an expression 
  75. # depending on q,qd,qdd 
  76. # used to compute Euler equations 
  77. #---------------------------------------------------------------
  78.  
  79.     ttvar:=proc(xx)
  80.             if type(xx,`indexed`) 
  81.         then cat(op(0,xx),`d`)[op(xx)]
  82.         else cat(xx,`d`) fi
  83.     end:
  84.  
  85.  
  86. time_diff:=proc(phi,q,qd,qdd)
  87.     local phi_copy,k,ttvar:
  88.     # subtitution to specify that q,qd ,qdd depends on time 
  89.     phi_copy:=phi:
  90.     phi_copy:=subs(map( xx-> xx=xx(t),[op(q),op(qd)]),phi_copy):
  91.     diff_phi:=diff(phi_copy,t):
  92.     # subtitution to come back to our variables 
  93.     diff_phi:=subs(map(xx->diff(xx(t),t)=ttvar(xx),[op(q),op(qd)]),
  94.             diff_phi):
  95.     diff_phi:=subs(map(xx->xx(t)=xx,[op(q),op(qd),op(qdd)]),diff_phi):
  96. end:
  97.  
  98.  
  99. #-----------------------------------------------------
  100. # Rewritting the Euler equations to have a canonical form 
  101. #           ..         .
  102. # El= ME(q)  q + CE(q) q^2 + RE(q)
  103. # Computation of ME,CE,RE 
  104. # CEuler returns a list [ME,CE,RE];
  105. #-----------------------------------------------------
  106.  
  107. CEuler:=proc(E,q,qd,qdd)
  108.     local Me,Ce,Re:
  109.     Me:=MME(E,q,qd,qdd):
  110.     Ce:=CCE(E,q,qd,qdd):
  111.     Re:=RRE(E,Me,Ce,q,qd,qdd):
  112.     [eval(Me),eval(Ce),eval(Re)]:
  113.     end:
  114.  
  115. MME:=proc(E,q,qd,qdd)
  116.     local E1:
  117.     E1:=eval(E):
  118.     genmatrix([seq(E1[i,1],i=1..nops(qdd))],qdd):
  119.     end:
  120.  
  121. #-----------------------------------------------------
  122. # extract the CE(q) matrix  El= ME(q)  qdd + CE(q) qd^2 + RE(q)
  123. #-----------------------------------------------------
  124.     ttvarp:=proc(xx)
  125.             if type(xx,`indexed`) 
  126.         then cat(op(0,xx),`2`)[op(xx)]
  127.         else cat(xx,`2`) fi
  128.     end:
  129.  
  130. CCE:=proc(E,q,qd,qdd)
  131.      local E_copy,q2d:
  132.     E_copy:=eval(E):
  133.     q2d:= map(x-> ttvarp(x),qd): 
  134.     E_copy:=subs( map(x-> x**2=ttvarp(x),qd),eval(E_copy));
  135.     genmatrix([seq(E_copy[i,1],i=1..nops(qd))],q2d);
  136.     end:
  137.  
  138. #-----------------------------------------------------
  139. # extract the RE(q) matrix  El= ME(q)  qdd + CE(q) qd^2 + RE(q)
  140. #-----------------------------------------------------
  141.  
  142. RRE:=proc(E,ME,CE,q,qd,qdd)
  143.     local MM:
  144.     MM:=add(    E,
  145.         add(multiply(ME,matrix(nops(q),1,qdd)),
  146.              multiply(CE,matrix(nops(q),1,map(x->x**2,qd)))),
  147.         1,-1);
  148.     MM:=map((x)-> simplify(x),eval(MM)):
  149.     end:
  150.  
  151. #
  152. #-----------------------------------------------------
  153. # FORTRAN GENERATION 
  154. #-----------------------------------------------------
  155. read ``.libname.`/macrofort.m`;
  156.  
  157. Gener:=proc(filename,fortranlist)
  158.     init_genfor();
  159.     optimized:=true:
  160.     writeto(filename):
  161.     genfor(flist):
  162.     writeto(terminal):
  163.     end:
  164.  
  165.